##Tugba Bozcaga and Kelly Pasolli
##Gov 2001 Project Graphs
my.d$lostXincome<- my.d$lostjob * my.d$income2
my.d$lessXincome<- my.d$HH_muchlessincome * my.d$income2
my.d$secureXincome<- my.d$lessecure * my.d$income2
my.d$male<-ifelse(my.d$gender_=="male", 1, 0)
my.d$domestic<-ifelse(my.d$marriage2=="Domestic partnership", 1, 0)
my.d$married<-ifelse(my.d$marriage2=="Married", 1, 0)
my.d$seperated<-ifelse(my.d$marriage2=="Seperated", 1, 0)
my.d$single<-ifelse(my.d$marriage2=="Single", 1, 0)
my.d$widowed<-ifelse(my.d$marriage2=="Widowed", 1, 0)
my.d$black<-ifelse(my.d$race2=="Black", 1, 0)
my.d$hispanic<-ifelse(my.d$race2=="Hispanic", 1, 0)
my.d$middle<-ifelse(my.d$race2=="Middle Eastern", 1, 0)
my.d$mixed<-ifelse(my.d$race2=="Mixed", 1, 0)
my.d$native<-ifelse(my.d$race2=="Native American", 1, 0)
my.d$other<-ifelse(my.d$race2=="Other", 1, 0)
my.d$white<-ifelse(my.d$race2=="White", 1, 0)
write.dta(my.d,file="data.final.dta")
c<-c(28, 29, 9, 11, 13, 14, 27, 4, 16, 15, 17, 21, 26, 7, 33, 34, 35, 36, 37, 40, 38,39, 41, 42,
     43, 44, 45, 12,25, 30,31,32)
c.mat<-my.d[,c]
##################################################################################################
###########CALCULATING EXPECTED VALUES OF Y FOR DIFFERENT AGE AND INCOME LEVELS###################
###########INCOME
library(Zelig)
z.out5i <- zelig(bi_welfare ~ 
                   dems + 
                   reps + 
                   lostjob +
                   HH_muchlessincome +
                   lessecure +
                   splostjob_ + income2 +
                   prevwelfa_ + stillunemployed + newjob + 
                   notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                   race2 + HHincnotreported + wave2 +
                   lostjob : income2 + HH_muchlessincome : income2 + lessecure : income2  
                   ,model = "probit", data=my.d, id = "caseid", robust = TRUE)
coefs<-summary(z.out5i)$coefficients[,1]
sigma<-summary(z.out5i)$cov.scaled

set.seed(1234)
library(mvtnorm)
sim.betas <- rmvnorm(n = 10000,
                   mean = coefs,
                   sigma = sigma)
Xc<-c()
for(i in 1:32){
  column.without.na<-na.omit(c.mat[,i])
  Xc[i]<-mean(column.without.na)
}

presvals <- seq(3,6,0.1)
ev.ests <- matrix(data = NA, ncol = length(presvals),
                  nrow=10000)

for(j in 1:length(presvals)){
  Xc.new <- Xc
  Xc.new[8] <- presvals[j]
  Xc.new[31]<- (presvals[j]* 0.042933032)
  Xc.new[32]<- (presvals[j]* 0.145477179)
  Xc.new[33]<- (presvals[j]* 0.239536016)
  ev.ests[,j] <- pnorm(Xc.new%*%t(sim.betas))
}

plot(presvals, apply(ev.ests,2,mean), ylim=c(0.1,0.5),xlab="Household Income(log)", 
     ylab="Expected Values",main="Expected Probability of Welfare Support by Income", cex.main=0.8)
segments(x0 = presvals, x1 = presvals,
         y0 = apply(ev.ests, 2, quantile, .025),
         y1 = apply(ev.ests, 2, quantile, .975))

##################################################################################################
###########CALCULATING EXPECTED VALUES OF Y FOR DIFFERENT AGE AND INCOME LEVELS###################
library(Zelig)
z.out5y <- zelig(bi_welfare ~ 
                   dems + 
                   reps + 
                   lostjob +
                   HH_muchlessincome +
                   lessecure +
                   splostjob_ +
                   prevwelfa_ + stillunemployed + newjob + 
                   notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
                   race2 + HHincnotreported + wave2 +
                   lostjob : age + HH_muchlessincome : age + lessecure : age  
                 ,model = "probit", data=my.d, id = "caseid", robust = TRUE)
coefs<-summary(z.out5y)$coefficients[,1]
sigma<-summary(z.out5y)$cov.scaled

set.seed(1234)
library(mvtnorm)
sim.betas <- rmvnorm(n = 10000,
                     mean = coefs,
                     sigma = sigma)
Xc<-c()
for(i in 1:32){
  column.without.na<-na.omit(c.mat[,i])
  Xc[i]<-mean(column.without.na)
}

presvals <- seq(15,95,1)
ev.ests <- matrix(data = NA, ncol = length(presvals),
                  nrow=10000)
Xc<-as.vector(c(1,Xc[-7]))
for(j in 1:length(presvals)){
  Xc.new <- Xc
  Xc.new[14] <- presvals[j]
  Xc.new[30]<- (presvals[j]* 0.042933032)
  Xc.new[31]<- (presvals[j]* 0.145477179)
  Xc.new[32]<- (presvals[j]* 0.239536016)
  ev.ests[,j] <- pnorm(Xc.new%*%t(sim.betas))
}

plot(presvals, apply(ev.ests,2,mean), ylim=c(0.2,0.8),xlab="Age", 
     ylab="Expected Values",main="Expected Probability of Welfare Support by Age", cex.main=0.8)
segments(x0 = presvals, x1 = presvals,
         y0 = apply(ev.ests, 2, quantile, .025),
         y1 = apply(ev.ests, 2, quantile, .975))






####################### GRAPHS #######################


#### Income and Job Loss ####

glm1<- glm(bi_welfare ~ 
             dems + 
             reps + 
             lostjob * income2 + 
             prevwelfa_ + stillunemployed + newjob + 
             notinlabormarket + logusd + ed2 + age + gender_ + marriage2 + 
             race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
           control  = list()) 

glm1
cov.mat<-vcov(glm1)
cov.mat
mod_frame=model.frame(glm1)
beta1<--.167
beta3<-.137
income.vals<-c(seq(2,6,0.01))
effects = beta1 + beta3*income.vals
var<-cov.mat[4,4] + (income.vals^2)*cov.mat[28,28] + 2*(income.vals)*cov.mat[4,28]
se<-sqrt(var)
z_score = qnorm(1 - ((1 - .95)/2))
upper=effects+z_score*se
lower=effects-z_score*se


plot(income.vals, effects, type="l", main="Figure 2: Relationship between Job Loss and
     Income on Support for Social Policy", xlab="Logged Household Income", ylab="Marginal Effect of Job Loss",
     cex.main=0.8, cex.lab=0.7, ylim=c(0, 0.6))
lines(y=upper, x=income.vals, lty=2)
lines(y=lower, x=income.vals, lty=2)

legend("topright",
       c("95% confidence intervals"),
       lty=c(1, 4), lwd=2, ncol=1, cex=0.5)

###### Income and Less Secure

glm2<- glm(bi_welfare ~ 
             dems + 
             reps +
             HH_muchlessincome * income2 + 
             prevwelfa_ + stillunemployed + newjob + notinlabormarket +
             logusd + ed2 + age + gender_ + marriage2 + 
             race2 + HHincnotreported + wave2, data=my.d, atmean = TRUE, robust = TRUE,clustervar1 =   "caseid",
           control  = list()) 

glm2
cov.mat<-vcov(glm2)
cov.mat
mod_frame=model.frame(glm2)
beta1<-.521
beta3<--.146
income.vals<-c(seq(2,6,0.01))
effects = beta1 + beta3*income.vals
var<-cov.mat[5,5] + (income.vals^2)*cov.mat[28,28] + 2*(income.vals)*cov.mat[5,28]
se<-sqrt(var)
z_score = qnorm(1 - ((1 - .95)/2))
upper=effects+z_score*se
lower=effects-z_score*se


plot(income.vals, effects, type="l", ylim=c(-0.7, 0.7), main="Figure 2: Relationship between Feeling Less
     Secure and Income on Support for Social Policy", xlab="Logged Household Income", ylab="Marginal Effect of Feeling Less Secure",
     cex.main=0.8, cex.lab=0.7)
lines(y=upper, x=income.vals, lty=2)
lines(y=lower, x=income.vals, lty=2)

legend("topright",
       c("95% confidence intervals"),
       lty=c(1, 4), lwd=2, ncol=1, cex=0.5)
